function  [Error_i, CS_i, W_i, W_experience_i,  Error_avg, CS_avg, W_avg, W_experience_avg, ext_i, ext_avg, ext_experience_i, ext_experience_avg, P_decision_j, P_experience_j,gvt_i,gvt_avg,SW_i,SW_avg,gvt_experience_i,gvt_experience_avg,SW_experience_i,SW_experience_avg]=Welfare_Heterweight_WTP_EEgap_parfor(param,param_FE,param_FEdemo,id_ij,J,p_elec,t_pigou,d_ext,elec,lifetime,r,MCPF,MT,price, price_denom, rebate_estar_denom, estar_denom, kwh, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom, LLC_5_18)

% Welfare formula:
% Delta CS = 1/eta *([log (sum_j U_j) +sum_j P_j(U^e-U)]_post_policy - [log (sum_j U_j) +sum_j P_j(U^e-U)]_pre_policy)

 %Parameters
  alpha=[0;param_FE(1:J-1,1)];
  eta=param(1);
  eta_experience=param(2);
  theta=param(3);
  tau=param(4);
  pi=param(5);
  beta=param_FEdemo;
  discount_fac=1/(1+r); 
  LLC=discount_fac*(1-discount_fac^lifetime)/(1-discount_fac); 
 
  FEprod_denom=alpha(prod_denom); 
 
 % kwh  
    kwh_denom=kwh(id_ij);
  
 % elec
    p_elec_ext=elec+p_elec+t_pigou;
  	elec_J=repmat(p_elec_ext,1,J);
    elec_J=elec_J';
    elec_denom=elec_J(id_ij).*kwh_denom;
    
       
 %Computations  
  				
 U_decision = FEprod_denom + eta_denom.*price_denom + pi*eta_denom.*rebate_estar_denom + tau*eta_denom.*estar_denom...
                + LLC_5_18*eta_denom.*theta_denom.*elec_denom...
  				+ beta(1)*type_id2_3_demo1_denom +beta(2)*type_id2_3_demo2_denom + beta(3)*type_id2_3_demo3_denom...
  				+ beta(4)*type_id2_3_demo4_denom +beta(5)*type_id2_3_demo5_denom + beta(6)*type_id2_3_demo6_denom...
				+ beta(7)*AV_demo1_denom +beta(8)*AV_demo2_denom + beta(9)*AV_demo3_denom...
  				+ beta(10)*AV_demo4_denom +beta(11)*AV_demo5_denom + beta(12)*AV_demo6_denom;				
 
 expU_decision=exp(U_decision);
  				
				
 U_experience= FEprod_denom + eta_experience_denom.*price_denom + LLC_5_18*eta_experience_denom.*elec_denom...
 				+ tau*eta_experience_denom.*estar_denom + pi*eta_experience_denom.*rebate_estar_denom... 				
 				+ beta(1)*type_id2_3_demo1_denom +beta(2)*type_id2_3_demo2_denom + beta(3)*type_id2_3_demo3_denom...
  				+ beta(4)*type_id2_3_demo4_denom +beta(5)*type_id2_3_demo5_denom + beta(6)*type_id2_3_demo6_denom...
				+ beta(7)*AV_demo1_denom +beta(8)*AV_demo2_denom + beta(9)*AV_demo3_denom...
  				+ beta(10)*AV_demo4_denom +beta(11)*AV_demo5_denom + beta(12)*AV_demo6_denom;				
  					
  expU_experience=exp(U_experience);
    
  U_decision_ij = sparse(id_dense_i,id_dense_j,U_decision,id_dense_m,id_dense_n)';
  num_decision = sparse(id_dense_i,id_dense_j,expU_decision,id_dense_m,id_dense_n)';
  denom_decision=sum(num_decision,2);

  U_experience_ij = sparse(id_dense_i,id_dense_j,U_experience,id_dense_m,id_dense_n)';
  num_experience = sparse(id_dense_i,id_dense_j,expU_experience,id_dense_m,id_dense_n)';
  denom_experience=sum(num_experience,2);
  
  P_decision_NK_J=num_decision./repmat(denom_decision,1,J);
  P_experience_NK_J=num_experience./repmat(denom_experience,1,J);
  
  P_decision_weigth_NK_J=weight.*P_decision_NK_J;
  P_experience_weigth_NK_J=weight.*P_experience_NK_J;
  
  P_decision_j=sum(P_decision,1)/MT;
  P_experience_j=sum(P_experience,1)/MT;
  
  Error_i=1000*abs(1/eta_experience)*sum(P_decision.*(U_experience_ij - U_decision_ij),2);
  CS_i=1000*abs(1/eta_experience)*(log(denom_decision));
  W_i=CS_i+Error_i;
  W_experience_i=1000*abs(1/eta_experience)*(log(denom_experience));
  
  Error_avg=sum(Error_i)/MT;
  CS_avg=sum(CS_i)/MT;
  W_avg=sum(W_i)/MT;
  W_experience_avg=sum(W_experience_i)/MT;
   
  ext_ij=LLC*kwh*d_ext*1000;
  ext_i=sum(P_decision.*ext_ij',2);
  ext_avg=sum(ext_i)/MT;
  ext_experience_i=sum(P_experience.*ext_ij',2);
  ext_experience_avg=sum(ext_experience_i)/MT;
  
  gvt_ij=LLC*kwh*t_pigou*1000;
  gvt_i=sum(P_decision.*gvt_ij',2);
  gvt_avg=sum(gvt_i)/MT;
  gvt_experience_i=sum(P_experience.*gvt_ij',2);
  gvt_experience_avg=sum(gvt_experience_i)/MT;
  
  SW_i=W_i-ext_i+MCPF*gvt_i;
  SW_avg=sum(SW_i)/MT;
  
  SW_experience_i=W_experience_i-ext_experience_i+MCPF*gvt_experience_i;
  SW_experience_avg=sum(SW_experience_i)/MT;
end 











